#####################################
#### Assign necessary libraries. ####

starttime <- Sys.time()

library(survival)
library(plyr)
library(optmatch)
library(rcbalance)
library(rcbsubset)
library(car)
library(data.table)
library(splitstackshape)
library(RCurl)
library(MASS)
library(pROC)


source("/Users/sharpe/Desktop/Young Surgeons/code/apply_penalty.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/subtabfun.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/meantab_both.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/baltab_ys.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/smahal_func.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/penalize.near.exact.v2.R")

###############################################
#### Read in the match and risk datasets. ####


ortho_risk_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/ortho_risk_patients.csv")
risk_score_patients_trad <- subset(ortho_risk_patients, era_trad==1)
risk_score_patients_mod <- subset(ortho_risk_patients, era_mod==1)

ortho_match_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/ortho_match_patients.csv")
trad_match_patients <- subset(ortho_match_patients, era_trad==1)
mod_match_patients <- subset(ortho_match_patients, era_mod==1)

NEW_PATIENTS_BEFORE <- subset(ortho_match_patients, new_surgeon_flag==1)
exp_patients_before <- subset(ortho_match_patients, new_surgeon_flag==0)

###############################################
###############################################

ortho_npgs <- c("hip_repair", "hip_replacement", "hip_replacement_revision", "knee_replacement", "knee_replacement_revision")
ortho_npgs_for_prop_score <- c("hip_repair", "hip_replacement", "hip_replacement_revision", "knee_replacement")


comorbidities <- c('silb_AbdCancer',  'silb_ANGINA', 'silb_ASTHMA', 'silb_CANCER', 'silb_CHF', 'silb_CHRNLUNG', 'silb_ChrPepticUlcer',
                   'silb_COAG', 'silb_CollagenVascular', 'silb_CongenitalCoag', 'silb_CUSHINGS', 'silb_DEMENT', 'silb_DIABETES', 'silb_GRAVES', 'silb_HTN', 'silb_HYPOTHY', 'silb_LIVER',
                   'silb_PARAPLEG', 'silb_PastARR', 'silb_PastMi', 'silb_PostPulmFibr', 'silb_RENLDYS', 'silb_renlfail',
                   'silb_SEIZUR', 'silb_SickleCell', 'silb_STROKE', 'silb_THROMB', 'silb_UnstAngina', 'silb_Valvulardis', 'silb_AIDS')
aids_position <- match("silb_AIDS", comorbidities)
#### Remove AIDS from the comorbidities, since it is pretty rare.
comorbidities_for_prop_score <- comorbidities[-aids_position]

prop_score_pat_vars_trad <- c("new_surgeon_flag", "admsn_in_quarters")
prop_score_pat_vars_mod <- c("new_surgeon_flag","admsn_in_quarters")



prop_score_pat_trad <- glm(formula = new_surgeon_flag ~ ., data=trad_match_patients[,c(prop_score_pat_vars_trad)], family=binomial)


prop_score_pat_mod_post_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==1),c(prop_score_pat_vars_mod)], family=binomial)
prop_score_pat_mod_pre_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==0),c(prop_score_pat_vars_mod)], family=binomial)

summary(prop_score_pat_trad)
summary(prop_score_pat_mod_post_2010)
summary(prop_score_pat_mod_pre_2010)


trad_match_patients$prop_prob_pat_trad <- predict(object=prop_score_pat_trad, newdata=trad_match_patients, type="response")
trad_match_patients$prop_log_pat_trad <- predict(object=prop_score_pat_trad, newdata=trad_match_patients, type="link", link="logit")

mod_match_patients$prop_prob_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"),0)

mod_match_patients$prop_prob_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"),0)

#### Calculate the standard deviation for the propensity scores. ####

prop_score_sd_trad <- sd(trad_match_patients[, c("prop_log_pat_trad")])
prop_score_sd_mod_post_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==1), c("prop_log_pat_mod_post_2010")])
prop_score_sd_mod_pre_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==0), c("prop_log_pat_mod_pre_2010")])

###############################################################################
#### Create a variable to prevent matrix/dataframe from becoming a vector. ####

trad_match_patients$for_matrix <- 1
mod_match_patients$for_matrix <- 1



#############################
#### Match the patients. ####


trad_match_patients <- trad_match_patients[order(-c(trad_match_patients$treatment)),]
mod_match_patients <- mod_match_patients[order(-c(mod_match_patients$treatment)),]

trad_treat <- trad_match_patients$treatment
mod_treat <- mod_match_patients$treatment


maha_vars_pat_unimp <- c("admsn_in_quarters", "for_matrix")
maha_vars_pat_imp <- c("admsn_in_quarters", "for_matrix")

trad_match_patients$exact_group_2 <- paste(trad_match_patients$pair_id_surg, sep="_")
mod_match_patients$exact_group_2 <- paste(mod_match_patients$pair_id_surg, mod_match_patients$oct_2010_plus, sep="_")



exact_vars_pat <- c("exact_group_2")
caliper_var_pat_trad <- c("admsn_in_quarters", "prop_log_pat_trad")
caliper_var_pat_mod <- c("admsn_in_quarters", "for_matrix", "prop_log_pat_mod_pre_2010", "prop_log_pat_mod_post_2010")

near_exact_vars <- NA



mnf.names <- vector("list",length=1)
mnf.names[[1]] <- c("admit_year")


maha_data_pat_unimp_trad <- trad_match_patients[,maha_vars_pat_unimp]
maha_data_pat_unimp_mod <- mod_match_patients[,maha_vars_pat_unimp]

maha_data_pat_imp_trad <- trad_match_patients[,maha_vars_pat_imp]
maha_data_pat_imp_mod <- mod_match_patients[,maha_vars_pat_imp]
exact_matrix_pat_trad <- trad_match_patients[,exact_vars_pat]
exact_matrix_pat_mod <- mod_match_patients[,exact_vars_pat]
calip_data_pat_trad <- trad_match_patients[,caliper_var_pat_trad]
calip_data_pat_mod <- mod_match_patients[,caliper_var_pat_mod]



################################################################
#### Create the important distance matrices; add penalties. ####

dist_struc_pats_imp_trad <- build.dist.struct(z=trad_treat, X = maha_data_pat_imp_trad, exact=exact_matrix_pat_trad, calip.option="none")

avg_imp_trad <- mean(unlist(dist_struc_pats_imp_trad))

for(i in 1:length(dist_struc_pats_imp_trad)){
  for(j in 1:length(dist_struc_pats_imp_trad[[i]])){
    dist_struc_pats_imp_trad[[i]][j] <- (dist_struc_pats_imp_trad[[i]][j]/avg_imp_trad)*1500
  }
  names(dist_struc_pats_imp_trad[[i]]) <- names(dist_struc_pats_imp_trad[[i]])
}

max(unlist(dist_struc_pats_imp_trad))

####

dist_struc_pats_imp_mod <- build.dist.struct(z=mod_treat, X = maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, calip.option="none")

avg_imp_mod <- mean(unlist(dist_struc_pats_imp_mod))

for(i in 1:length(dist_struc_pats_imp_mod)){
  for(j in 1:length(dist_struc_pats_imp_mod[[i]])){
    dist_struc_pats_imp_mod[[i]][j] <- (dist_struc_pats_imp_mod[[i]][j]/avg_imp_mod)*1500
  }
  names(dist_struc_pats_imp_mod[[i]]) <- names(dist_struc_pats_imp_mod[[i]])
}

max(unlist(dist_struc_pats_imp_mod))


####

dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=maha_data_pat_imp_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_trad, cal.size=0.20, pen.size=2000, var="prop_log_pat_trad")

####

dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_post_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_pre_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_pre_2010")

####



dist_struc_pats_trad <- dist_struc_pats_imp_trad
dist_struc_pats_mod <- dist_struc_pats_imp_mod

####

treated_info_pats_trad <- subset(trad_match_patients, treatment==1)
control_info_pats_trad <- subset(trad_match_patients, treatment==0)

treated_info_pats_mod <- subset(mod_match_patients, treatment==1)
control_info_pats_mod <- subset(mod_match_patients, treatment==0)


patient_match_trad <- rcbsubset(distance.structure = dist_struc_pats_trad, fb.list=mnf.names, treated.info = treated_info_pats_trad, control.info = control_info_pats_trad,
                                penalty=2, tol=10)

patient_match_mod <- rcbsubset(distance.structure = dist_struc_pats_mod, fb.list=mnf.names,treated.info = treated_info_pats_mod, control.info = control_info_pats_mod,
                               penalty=2, tol=10)


tem_match_pats_trad <- treated_info_pats_trad[as.numeric(rownames(patient_match_trad$matches)),]
sda_match_pats_trad <- control_info_pats_trad[(patient_match_trad$matches),]
pair_id_pats_trad <- c(seq(1:dim(tem_match_pats_trad)[1]), seq(1:dim(sda_match_pats_trad)[1]))
combined_pats_trad <- rbind(tem_match_pats_trad, sda_match_pats_trad)
combined_pats_trad$pair_id_pats <- pair_id_pats_trad

tem_match_pats_mod <- treated_info_pats_mod[as.numeric(rownames(patient_match_mod$matches)),]
sda_match_pats_mod <- control_info_pats_mod[(patient_match_mod$matches),]
pair_id_pats_mod <- c(seq(1:dim(tem_match_pats_mod)[1]), seq(1:dim(sda_match_pats_mod)[1]))
combined_pats_mod <- rbind(tem_match_pats_mod, sda_match_pats_mod)
combined_pats_mod$pair_id_pats <- pair_id_pats_mod

#### Add propensity scores to the before datasets.

NEW_PATIENTS_BEFORE$prop_prob_pat_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=prop_score_pat_trad, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=prop_score_pat_trad, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)



exp_patients_before$prop_prob_pat_trad <- ifelse(exp_patients_before$era_trad==1, predict(object=prop_score_pat_trad, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_trad <- ifelse(exp_patients_before$era_trad==1, predict(object=prop_score_pat_trad, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)

NEW_PATIENTS_BEFORE$treatment <- 1



###############################################
#### Save necessary pieces of information. ####

#### Save the before match datasets. ####
write.csv(x=NEW_PATIENTS_BEFORE, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_new_pats.csv")
write.csv(x=exp_patients_before, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_exp_pats.csv")



#### Save the matched data. ####
write.csv(combined_pats_trad,  "/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_matched_trad_pats.csv")
write.csv(combined_pats_mod,  "/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_matched_mod_pats.csv")

#### Save the models. ####
prop_score_trad_summary <- summary(prop_score_pat_trad)
write.csv(prop_score_trad_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_propensity_model_trad_ortho.csv")

prop_score_mod_summary <- summary(prop_score_pat_mod_pre_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_propensity_model_mod_ortho_pre_2010.csv")
prop_score_mod_summary <- summary(prop_score_pat_mod_post_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_propensity_model_mod_ortho_post_2010.csv")



######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################

rm(list=ls())

######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################
######################################################################################################################################################

#####################################
#### Assign necessary libraries. ####

starttime <- Sys.time()

library(survival)
library(plyr)
library(optmatch)
library(rcbalance)
library(rcbsubset)
library(car)
library(data.table)
library(splitstackshape)
library(RCurl)
library(MASS)
library(pROC)


source("/Users/sharpe/Desktop/Young Surgeons/code/apply_penalty.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/subtabfun.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/meantab_both.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/baltab_ys.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/smahal_func.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/penalize.near.exact.v2.R")

###############################################
#### Read in the match and risk datasets. ####


gensurg_risk_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/gensurg_risk_patients.csv")
risk_score_patients_trad <- subset(gensurg_risk_patients, era_trad==1)
risk_score_patients_mod <- subset(gensurg_risk_patients, era_mod==1)

gensurg_match_patients <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/gensurg_match_patients.csv")
trad_match_patients <- subset(gensurg_match_patients, era_trad==1)
mod_match_patients <- subset(gensurg_match_patients, era_mod==1)

NEW_PATIENTS_BEFORE <- subset(gensurg_match_patients, new_surgeon_flag==1)
exp_patients_before <- subset(gensurg_match_patients, new_surgeon_flag==0)


################################
#### Define some variables. ####


gensurg_npgs <- c("cholecystectomy", "hernia_abdominal_open", "colectomy_partial__open", "appendectomy",                               
                  "lysis_of_adhesions", "mastectomy", "proctopexy",
                  "hernia_groin_open", "small_bowel_other", "pyloroplasty",
                  "proctectomy", "stomach_antireflux", "large_bowel_other",
                  "thyroidectomy_partial", "biliary_other", "small_bowel_resection",
                  "pd_access_procedure", "ulcer", "pancreatectomy",
                  "parathyroidectomy", "esophagectomy", "adrenalectomy",
                  "stomach_gastric_bypass_nonbariatric", "stomach_other", "splenectomy",
                  "liver_partial_hepatectomy", "colectomy_total__open", "stomach_partial_gastrectomy",
                  "thyroidectomy_total", "stomach_total_gastrectomy", "bariatric",
                  "biliary_common_duct", "liver_other", "thyroidectomy_substernal",
                  "esophagomyotomy", "hernia_abdominal_lap", "colectomy_partial__lap",
                  "colectomy_total__lap", "hernia_groin_lap")
gensurg_npgs_for_prop_score <- gensurg_npgs[gensurg_npgs != "colectomy_partial__open"]

comorbidities <- c('silb_AbdCancer',  'silb_ANGINA', 'silb_ASTHMA', 'silb_CANCER', 'silb_CHF', 'silb_CHRNLUNG', 'silb_ChrPepticUlcer',
                   'silb_COAG', 'silb_CollagenVascular', 'silb_CongenitalCoag', 'silb_CUSHINGS', 'silb_DEMENT', 'silb_DIABETES', 'silb_GRAVES', 'silb_HTN', 'silb_HYPOTHY', 'silb_LIVER',
                   'silb_PARAPLEG', 'silb_PastARR', 'silb_PastMi', 'silb_PostPulmFibr', 'silb_RENLDYS', 'silb_renlfail',
                   'silb_SEIZUR', 'silb_SickleCell', 'silb_STROKE', 'silb_THROMB', 'silb_UnstAngina', 'silb_Valvulardis', 'silb_AIDS')
aids_position <- match("silb_AIDS", comorbidities)
#### Remove AIDS from the comorbidities, since it is pretty rare.
comorbidities_for_prop_score <- comorbidities[-aids_position]




prop_score_pat_vars_trad <- c("new_surgeon_flag", "admsn_in_quarters")
prop_score_pat_vars_mod <- c("new_surgeon_flag","admsn_in_quarters")



prop_score_pat_trad <- glm(formula = new_surgeon_flag ~ ., data=trad_match_patients[,c(prop_score_pat_vars_trad)], family=binomial)


prop_score_pat_mod_post_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==1),c(prop_score_pat_vars_mod)], family=binomial)
prop_score_pat_mod_pre_2010 <- glm(formula = new_surgeon_flag ~ ., data=mod_match_patients[which(mod_match_patients$oct_2010_plus==0),c(prop_score_pat_vars_mod)], family=binomial)

summary(prop_score_pat_trad)
summary(prop_score_pat_mod_post_2010)
summary(prop_score_pat_mod_pre_2010)


trad_match_patients$prop_prob_pat_trad <- predict(object=prop_score_pat_trad, newdata=trad_match_patients, type="response")
trad_match_patients$prop_log_pat_trad <- predict(object=prop_score_pat_trad, newdata=trad_match_patients, type="link", link="logit")

mod_match_patients$prop_prob_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"),0)

mod_match_patients$prop_prob_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"),0)

#### Calculate the standard deviation for the propensity scores. ####

prop_score_sd_trad <- sd(trad_match_patients[, c("prop_log_pat_trad")])
prop_score_sd_mod_post_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==1), c("prop_log_pat_mod_post_2010")])
prop_score_sd_mod_pre_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==0), c("prop_log_pat_mod_pre_2010")])

###############################################################################
#### Create a variable to prevent matrix/dataframe from becoming a vector. ####

trad_match_patients$for_matrix <- 1
mod_match_patients$for_matrix <- 1




#############################
#### Match the patients. ####


trad_match_patients <- trad_match_patients[order(-c(trad_match_patients$treatment)),]
mod_match_patients <- mod_match_patients[order(-c(mod_match_patients$treatment)),]

trad_treat <- trad_match_patients$treatment
mod_treat <- mod_match_patients$treatment


maha_vars_pat_unimp <- c("admsn_in_quarters", "for_matrix")
maha_vars_pat_imp <- c("admsn_in_quarters", "for_matrix")

trad_match_patients$exact_group_2 <- paste(trad_match_patients$pair_id_surg, sep="_")
mod_match_patients$exact_group_2 <- paste(mod_match_patients$pair_id_surg, mod_match_patients$oct_2010_plus, sep="_")



exact_vars_pat <- c("exact_group_2")
caliper_var_pat_trad <- c("admsn_in_quarters", "prop_log_pat_trad")
caliper_var_pat_mod <- c("admsn_in_quarters", "for_matrix", "prop_log_pat_mod_pre_2010", "prop_log_pat_mod_post_2010")

near_exact_vars <- NA



# mnf.names <- vector("list",length=1)
# mnf.names[[1]] <- c("admit_year")


maha_data_pat_unimp_trad <- trad_match_patients[,maha_vars_pat_unimp]
maha_data_pat_unimp_mod <- mod_match_patients[,maha_vars_pat_unimp]

maha_data_pat_imp_trad <- trad_match_patients[,maha_vars_pat_imp]
maha_data_pat_imp_mod <- mod_match_patients[,maha_vars_pat_imp]
exact_matrix_pat_trad <- trad_match_patients[,exact_vars_pat]
exact_matrix_pat_mod <- mod_match_patients[,exact_vars_pat]
calip_data_pat_trad <- trad_match_patients[,caliper_var_pat_trad]
calip_data_pat_mod <- mod_match_patients[,caliper_var_pat_mod]



################################################################
#### Create the important distance matrices; add penalties. ####

dist_struc_pats_imp_trad <- build.dist.struct(z=trad_treat, X = maha_data_pat_imp_trad, exact=exact_matrix_pat_trad, calip.option="none")

avg_imp_trad <- mean(unlist(dist_struc_pats_imp_trad))

for(i in 1:length(dist_struc_pats_imp_trad)){
  for(j in 1:length(dist_struc_pats_imp_trad[[i]])){
    dist_struc_pats_imp_trad[[i]][j] <- (dist_struc_pats_imp_trad[[i]][j]/avg_imp_trad)*1500
  }
  names(dist_struc_pats_imp_trad[[i]]) <- names(dist_struc_pats_imp_trad[[i]])
}

max(unlist(dist_struc_pats_imp_trad))

####

dist_struc_pats_imp_mod <- build.dist.struct(z=mod_treat, X = maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, calip.option="none")

avg_imp_mod <- mean(unlist(dist_struc_pats_imp_mod))

for(i in 1:length(dist_struc_pats_imp_mod)){
  for(j in 1:length(dist_struc_pats_imp_mod[[i]])){
    dist_struc_pats_imp_mod[[i]][j] <- (dist_struc_pats_imp_mod[[i]][j]/avg_imp_mod)*1500
  }
  names(dist_struc_pats_imp_mod[[i]]) <- names(dist_struc_pats_imp_mod[[i]])
}

max(unlist(dist_struc_pats_imp_mod))


####

dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=maha_data_pat_imp_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_trad, cal.size=0.20, pen.size=2000, var="prop_log_pat_trad")

####

dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_post_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_pre_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_pre_2010")

####



dist_struc_pats_trad <- dist_struc_pats_imp_trad
dist_struc_pats_mod <- dist_struc_pats_imp_mod

####

treated_info_pats_trad <- subset(trad_match_patients, treatment==1)
control_info_pats_trad <- subset(trad_match_patients, treatment==0)

treated_info_pats_mod <- subset(mod_match_patients, treatment==1)
control_info_pats_mod <- subset(mod_match_patients, treatment==0)


patient_match_trad <- rcbsubset(distance.structure = dist_struc_pats_trad, fb.list=mnf.names, treated.info = treated_info_pats_trad, control.info = control_info_pats_trad,
                                penalty=2, tol=10)

patient_match_mod <- rcbsubset(distance.structure = dist_struc_pats_mod, fb.list=mnf.names, treated.info = treated_info_pats_mod, control.info = control_info_pats_mod,
                               penalty=2, tol=10)


tem_match_pats_trad <- treated_info_pats_trad[as.numeric(rownames(patient_match_trad$matches)),]
sda_match_pats_trad <- control_info_pats_trad[(patient_match_trad$matches),]
pair_id_pats_trad <- c(seq(1:dim(tem_match_pats_trad)[1]), seq(1:dim(sda_match_pats_trad)[1]))
combined_pats_trad <- rbind(tem_match_pats_trad, sda_match_pats_trad)
combined_pats_trad$pair_id_pats <- pair_id_pats_trad

tem_match_pats_mod <- treated_info_pats_mod[as.numeric(rownames(patient_match_mod$matches)),]
sda_match_pats_mod <- control_info_pats_mod[(patient_match_mod$matches),]
pair_id_pats_mod <- c(seq(1:dim(tem_match_pats_mod)[1]), seq(1:dim(sda_match_pats_mod)[1]))
combined_pats_mod <- rbind(tem_match_pats_mod, sda_match_pats_mod)
combined_pats_mod$pair_id_pats <- pair_id_pats_mod

#### Add propensity scores to the before datasets.

NEW_PATIENTS_BEFORE$prop_prob_pat_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=prop_score_pat_trad, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=prop_score_pat_trad, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)



exp_patients_before$prop_prob_pat_trad <- ifelse(exp_patients_before$era_trad==1, predict(object=prop_score_pat_trad, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_trad <- ifelse(exp_patients_before$era_trad==1, predict(object=prop_score_pat_trad, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)

NEW_PATIENTS_BEFORE$treatment <- 1




###############################################
#### Save necessary pieces of information. ####

#### Save the before match datasets. ####
write.csv(x=NEW_PATIENTS_BEFORE, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_new_pats.csv")
write.csv(x=exp_patients_before, file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_exp_pats.csv")



#### Save the matched data. ####
write.csv(combined_pats_trad,  "/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_matched_trad_pats.csv")
write.csv(combined_pats_mod,  "/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_matched_mod_pats.csv")

#### Save the models. ####
prop_score_trad_summary <- summary(prop_score_pat_trad)
write.csv(prop_score_trad_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_propensity_model_trad_gensurg.csv")

prop_score_mod_summary <- summary(prop_score_pat_mod_pre_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_propensity_model_mod_gensurg_pre_2010.csv")
prop_score_mod_summary <- summary(prop_score_pat_mod_post_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_propensity_model_mod_gensurg_post_2010.csv")





